﻿#region License Info
//Component of Cronos Package, http://www.codeplex.com/cronos
//Copyright (C) 2009 Anthony Brockwell

//This program is free software; you can redistribute it and/or
//modify it under the terms of the GNU General Public License
//as published by the Free Software Foundation; either version 2
//of the License, or (at your option) any later version.

//This program is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//GNU General Public License for more details.

//You should have received a copy of the GNU General Public License
//along with this program; if not, write to the Free Software
//Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#endregion


using System;
using System.Collections.Generic;
using System.Drawing;
using System.Windows.Forms;
using ABMath.Miscellaneous;
using ABMath.ModelFramework.Data;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.LinearAlgebra;

namespace ExtraControlLibrary
{
    public partial class MultiScaleCorrelogram : UserControl
    {
        private TimeSeries theTimeSeries;
        public TimeSeries TheTimeSeries
        {
            get { return theTimeSeries; } 
            set
            {
                theTimeSeries = value;
                RebuildAll();
            }
        }

        public int NumIterations
        { get; set; }

        private SurfaceEstimate Surface
        {
            get
            {
                return shadedSurfaceBox1.Surface;
            }
            set
            {
                if (value == null)
                {
                    shadedSurfaceBox1.Surface = null;
                    return;
                }
                if (value.Dimension != 2)
                    throw new ApplicationException("Surface must be defined on a two-dimensional domain.");
                shadedSurfaceBox1.Surface = value;
                Recompute();
            }
        }

        private void Recompute()
        {
            if (shadedSurfaceBox1.Surface == null)
                return;
            
            shadedSurfaceBox1.Recompute();

            var legend = new Dictionary<Brush, double>();
            double min = shadedSurfaceBox1.MinValue;
            double max = shadedSurfaceBox1.MaxValue;
            int numDivs = 7;
            for (int i = 0; i < numDivs; ++i)
            {
                double value = i * (max - min) / (numDivs - 1) + min;
                legend.Add(shadedSurfaceBox1.GetBrushForValue(value), value);
            }

            shadedSurfaceBox1.Legend = legend;
        }

        private void RebuildAll()
        {
            if (theTimeSeries == null)
                return;

            int numSubits = 80;

            SurfaceEstimate correlations;

            double oneDayTicks = new TimeSpan(0, 24, 0, 0).Ticks;
            double oneMonthTicks = new TimeSpan(30, 0, 0, 0).Ticks;

            double oneMinuteTicks = new TimeSpan(0, 0, 1, 0).Ticks;
            double twoHoursTicks = new TimeSpan(0, 2, 0, 0).Ticks;
            double fortyDaysTicks = new TimeSpan(40, 0, 0, 0).Ticks;
            double t0 = theTimeSeries.TimeStamp(0).Ticks;
            double t1 = theTimeSeries.TimeStamp(theTimeSeries.Count - 1).Ticks;

            var lowerLimit = new Vector(new[] { Math.Log(oneMinuteTicks), Math.Log(oneMinuteTicks) });
            var upperLimit = new Vector(new[] { Math.Log(fortyDaysTicks), Math.Log(fortyDaysTicks) });

            correlations = new SurfaceEstimate(lowerLimit, upperLimit, 40);

            // now fill in autocorrelations by sampling
            var cud = new ContinuousUniformDistribution();

            var hs = new HaltonSequence(2);

            double totalGamma = 0;
            double totalGammaWeight = 0;

            var cache = theTimeSeries.BuildLookupCache(128);

            for (int iteration = 0; iteration < NumIterations; ++iteration)
            {
                // pick a random point
                Vector draw = hs.GetNext();
                double h = draw[0] * (upperLimit[0] - lowerLimit[0]) + lowerLimit[0];
                double b = draw[1] * (upperLimit[1] - lowerLimit[1]) + lowerLimit[1];

                double realtimeh = Math.Exp(h);
                double hindays = realtimeh / oneDayTicks;
                double realtimeb = Math.Exp(b);
                double bindays = realtimeb / oneDayTicks;

                double value = 0;
                for (int subits = 0; subits < numSubits; ++subits)
                {
                    double t = cud.NextDouble() * (t1 - t0 - realtimeh - realtimeb) + t0;

                    var times = new[] 
                    {   
                        new DateTime((long)t), 
                        new DateTime((long)(t+realtimeb)), 
                        new DateTime((long)(t+realtimeh)), 
                        new DateTime((long)(t+realtimeh + realtimeb)) 
                    };
                    var vals = theTimeSeries.ValuesAtTimes(times, cache);

                    double xtplushplusb = Math.Log(vals[3]);
                    double xtplush = Math.Log(vals[2]);
                    double xtplusb = Math.Log(vals[1]);
                    double xt = Math.Log(vals[0]);

                    if (h > b)
                        value += (xtplusb - xt) * (xtplushplusb - xtplush) / bindays;
                    else
                        value += (xtplushplusb - xtplusb) * (xtplusb - xt) / Math.Sqrt(bindays * hindays);

                    totalGamma += ((xtplusb - xt) * (xtplusb - xt)) / bindays;
                    totalGammaWeight += 1.0;
                }
                value /= numSubits;

                correlations.AddPoint(new Vector(new[] { h, b }), value, 1.0);
            }

            double gammahat = totalGamma / totalGammaWeight;
            correlations.MultiplyBy(1.0 / gammahat);
            Surface = correlations;

            // now save it all
            //using (var swriter = new StreamWriter(new FileStream("c:\\temp\\cgtest.out", FileMode.Create)))
            //{
            //    foreach (var x in testpoints)
            //    {
            //        foreach (var y in x)
            //            swriter.Write("{0:0.00000}\t", y);
            //        swriter.WriteLine();
            //    }
            //}
        }

        public MultiScaleCorrelogram()
        {
            TheTimeSeries = null;
            NumIterations = 500000;
            InitializeComponent();
        }
    }
}
